#include<math.h>
#include<GL/gl.h>
#include<GL/glu.h>
#include<GL/glut.h>
#include<iostream>
#include<vector>
#include<stdio.h>
#include<time.h>

using namespace std;

GLfloat screen_width,screen_height,aspect_ratio;
clock_t prev_time,next_time;

class Vector3D
{
public:
	GLfloat x,y,z;
	Vector3D(){
		x=y=z=0;
	}
	Vector3D(GLfloat X, GLfloat Y, GLfloat Z){
		x=X;
		y=Y;
		z=Z;
	}
	Vector3D operator+(Vector3D v){
		return Vector3D( x+v.x, y+v.y, z+v.z );
	}
	Vector3D operator-(Vector3D v){
		return Vector3D( x-v.x, y-v.y, z-v.z );
	}
	Vector3D operator*(GLfloat scalar){
		return Vector3D( x*scalar, y*scalar, z*scalar );
	}
	Vector3D operator/(GLfloat scalar){
		return Vector3D( x/scalar, y/scalar, z/scalar );
	}
	void operator=(Vector3D replacement)
	{
		x=replacement.x;
		y=replacement.y;
		z=replacement.z;
	}
	GLfloat get_magnitude(){
		return sqrt(pow(x,2)+pow(y,2)+pow(z,2));
	}
	void normalize(){
		Vector3D temp;
		temp = (*this) / get_magnitude();

		x=temp.x;
		y=temp.y;
		z=temp.z;

	}
	void rotate(GLfloat angle, Vector3D axis){

		GLfloat s=sin(angle);
		GLfloat c=cos(angle);

		axis.normalize();
		GLfloat prev_x = x, prev_y = y, prev_z = z;
		GLfloat u = axis.x, v = axis.y, w = axis.z;

		x = u *(u*prev_x + v*prev_y + w*prev_z) + c *(prev_x *(v*v + w*w) - u*(v*prev_y + w*prev_z)) + s*(-w*prev_y + v*prev_z);
		y = v *(u*prev_x + v*prev_y + w*prev_z) + c *(prev_y *(u*u + w*w) - v*(u*prev_x + w*prev_z)) + s*(w*prev_y - u*prev_z);
		z = w *(u*prev_x + v*prev_y + w*prev_z) + c *(prev_z *(u*u + v*v) - w*(u*prev_x + v*prev_y)) + s*(-v*prev_x + u*prev_y);
	}
};

class Point3D
{
public:
	GLfloat x,y,z;
	Point3D(){
		x=y=z=0;
	}
	Point3D(GLfloat X, GLfloat Y, GLfloat Z){
		x=X;
		y=Y;
		z=Z;
	}
	Point3D operator+(Vector3D p)
	{
		return Point3D(x+p.x, y+p.y, z+p.z);
	}
	Vector3D operator-(Point3D p)
	{
		return Vector3D(x-p.x, y-p.y, z-p.z);
	}
	Point3D operator=(Point3D replacement)
	{
		x=replacement.x;
		y=replacement.y;
		z=replacement.z;
	}
};

class Affine
{
private:
	GLfloat transformation_matrix[4][4];
	void multiply(Affine matrix)
	{
		GLfloat temp[4][4];

		int i,j,k;
		for( i=0; i<4; i++ )
			for( j=0; j<4; j++ )
				temp[i][j]=transformation_matrix[i][j];

		for( i=0; i<4; i++ )
			for( j=0; j<4; j++ )
				for( transformation_matrix[i][j]=0, k=0; k<4; k++ )
					transformation_matrix[i][j] += matrix.transformation_matrix[i][k] * temp[k][j];
	}
public:
	Affine()
	{
		set_identity();
	}
	void print_affine()
	{
		cout<<"\nstart\n";
		for(int i=0;i<4;i++){
			cout<<"\n";
			for(int j=0;j<4;j++)
				cout<<transformation_matrix[i][j]<<" ";
		}
		cout<<"\nend\n";
	}
	void set_identity()
	{
		int i,j;
		for( i=0; i<4; i++ )
			for( j=0; j<4; j++ )
				if( i==j )
					transformation_matrix[i][j]=1;
				else
					transformation_matrix[i][j]=0;
	}
	void apply_transformation()
	{
		glMultMatrixf(transformation_matrix[0]);
	}

	void translate(Vector3D displacement)
	{
		Affine temp;
		temp.transformation_matrix[3][0] = displacement.x;
		temp.transformation_matrix[3][1] = displacement.y;
		temp.transformation_matrix[3][2] = displacement.z;
		multiply(temp);
	}
	void rotate(GLfloat angle, Vector3D axis)
	{
		axis.normalize();
		Affine temp;
		GLfloat c = cos(angle);
		GLfloat s = -sin(angle);
		GLfloat mc = 1 - c;
		temp.transformation_matrix[0][0] = c + mc * axis.x * axis.x;
		temp.transformation_matrix[1][0] = mc * axis.x * axis.y + s * axis.z;
		temp.transformation_matrix[2][0] = mc * axis.x * axis.z - s * axis.y;
		temp.transformation_matrix[0][1] = mc * axis.y * axis.x - s * axis.z;
		temp.transformation_matrix[1][1] = c + mc * axis.y * axis.y;
		temp.transformation_matrix[2][1] = mc * axis.y * axis.z + s * axis.x;
		temp.transformation_matrix[0][2] = mc * axis.z * axis.x + s * axis.y;
		temp.transformation_matrix[1][2] = mc * axis.z * axis.y - s * axis.x;
		temp.transformation_matrix[2][2] = c + mc * axis.z * axis.z;
		multiply(temp);

	}
	void scale(GLfloat sx,GLfloat sy, GLfloat sz)
	{
		Affine temp;
		temp.transformation_matrix[0][0] = sx;
		temp.transformation_matrix[1][1] = sy;
		temp.transformation_matrix[2][2] = sz;
		multiply(temp);

	}
};